package edu.cornell.lassp.houle.RngPack;

import java.io.Serializable;
import java.util.Date;

//
// RngPack 1.1a by Paul Houle
// http://www.honeylocust.com/RngPack/
//

/**
 * 
 * <TT>RANMAR</TT> is a lagged Fibonacci generator proposed by Marsaglia and Zaman and is a good
 * research grade generator. This version of <TT>RANMAR</TT> is based on the paper by James, which
 * is a good reference for the properties of <TT>RANMAR</TT> and several other generators.
 * 
 * <BR>
 * <B>REFERENCES:</B> <BR>
 * F. James, <CITE>Comp. Phys. Comm.</CITE> <STRONG>60</STRONG> (1990) p 329-344 <BR>
 * and was originally described in <BR>
 * G. Marsaglia, A. Zaman and W.-W Tsang, <CITE>Stat. Prob. Lett</CITE> <STRONG>9</STRONG> (1990)
 * p 35.
 * 
 * 
 * <P>
 * <A HREF="/RngPack/src/edu/cornell/lassp/houle/RngPack/Ranmar.java"> Source code </A> is
 * available.
 * 
 * @author <A HREF="http://www.honeylocust.com/"> Paul Houle </A> (E-mail: <A
 *         HREF="mailto:paul@honeylocust.com">paul@honeylocust.com</A>)
 * @version 1.1a
 */

public class Ranmar
    extends RandomSeedable implements Serializable {

  double c, cd, cm, u[], uvec[];
  int i97, j97;

  /**
   * Default seed. <CODE>DEFSEED=54217137</CODE>
   */
  public static int DEFSEED = 54217137;

  /**
   * The 46,009,220nd prime number, he largest prime less than 9*10<SUP>8</SUP>. Used as a modulus
   * because this version of <TT>RANMAR</TT> needs a seed between 0 and 9*10<SUP>8</SUP> and
   * <CODE>BIG_PRIME</CODE> isn't commensurate with any regular period. <CODE>BIG_PRIME=
   * 899999963</CODE>
   */
  public static int BIG_PRIME = 899999963;

  /**
   * 
   * Initialize Ranmar with a specified integer seed
   * 
   * @param ijkl seed integer; <TT>Ranmar(int ijkl)</TT> takes uses <TT>ijkl</TT> modulus <TT>BIG_PRIME</TT>
   *        as a seed for <TT>RANMAR.</TT>
   * 
   */

  public Ranmar(int ijkl) {
    ranmarin(Math.abs(ijkl % BIG_PRIME));
  };

  /**
   * 
   * Initialize Ranmar with a specified long seed
   * 
   * @param ijkl seed long; <TT>Ranmar(long ijkl)</TT> takes uses <TT>ijkl</TT> modulus <TT>BIG_PRIME</TT>
   *        as a seed for <TT>RANMAR.</TT>
   * 
   */

  public Ranmar(long ijkl) {
    ranmarin((int) Math.abs(ijkl % BIG_PRIME));
  };

  /**
   * 
   * Initialize Ranmar with a default seed taken from Marsaglia and Zaman's paper. Equivalent to
   * <CODE>Ranmar(54217137).</CODE>
   * 
   */

  public Ranmar() {
    ranmarin(DEFSEED);
  };

  /**
   * 
   * Seed <TT>RANMAR</TT> from the clock.
   * 
   * <PRE>
   * 
   * RandomElement e=new Ranmar(new Date());
   * 
   * </PRE>
   * 
   * @param d a Date object to seed Ranmar with, typically <CODE>new Date()</CODE>
   * 
   */

  public Ranmar(Date d) {

    ranmarin((int) ClockSeed(d) % BIG_PRIME);
  };

  /**
   * 
   * Internal methods: ranmarin is the initialization code for the generator.
   * 
   */

  void ranmarin(int ijkl) {

    int ij, kl;
    int i, ii, j, jj, k, l, m;
    double s, t;

    u = new double[97];
    uvec = new double[97];

    ij = ijkl / 30082;
    kl = ijkl - 30082 * ij;

    i = ((ij / 177) % 177) + 2;
    j = (ij % 177) + 2;
    k = ((kl / 169) % 178) + 1;
    l = kl % 169;
    for (ii = 0; ii < 97; ii++) {
      s = 0.0;
      t = 0.5;
      for (jj = 0; jj < 24; jj++) {
        m = (((i * j) % 179) * k) % 179;
        i = j;
        j = k;
        k = m;
        l = (53 * l + 1) % 169;
        if (((l * m) % 64) >= 32)
          s += t;
        t *= 0.5;
      }
      u[ii] = s;
    }
    c = 362436.0 / 16777216.0;
    cd = 7654321.0 / 16777216.0;
    cm = 16777213.0 / 16777216.0;
    i97 = 96;
    j97 = 32;

  };

  /**
   * The generator
   * 
   * @return a pseudo random number
   */

  final public double raw() {

    double uni;

    uni = u[i97] - u[j97];
    if (uni < 0.0)
      uni += 1.0;
    u[i97] = uni;
    if (--i97 < 0)
      i97 = 96;
    if (--j97 < 0)
      j97 = 96;
    c -= cd;
    if (c < 0.0)
      c += cm;
    uni -= c;
    if (uni < 0.0)
      uni += 1.0;
    return (uni);
  };

  /**
   * 
   * A version of the generator for filling arrays, inlined for speed
   * 
   * @param d an array of doubles to be filled
   * @param n size of the array
   * 
   */

  final public void raw(double d[], int n) {

    double uni;

    for (int i = 0; i < n; i++) {
      uni = u[i97] - u[j97];
      if (uni < 0.0)
        uni += 1.0;
      u[i97] = uni;
      if (--i97 < 0)
        i97 = 96;
      if (--j97 < 0)
        j97 = 96;
      c -= cd;
      if (c < 0.0)
        c += cm;
      uni -= c;
      if (uni < 0.0)
        uni += 1.0;
      d[i] = uni;
    }
    ;
  };

};
